Journal of Computational Chemistry

Predicting the location of metal-binding sites in proteins is crucial for fundamental biological questions and biotechnological applications. Over the past decade, the rise in metal-bound protein structures in the Protein Data Bank, combined with advanced statistical models such as deep learning, has accelerated the development of metal-binding site prediction tools. Several approaches are now available, offering high-quality benchmarks and predictive performance. Our initial development in this area is BioMetAll, whose first version was based on backbone pre-organization. Here, we introduce its second version, featuring two major updates: 1) metal-specific scoring functions and 2) prediction using backbone geometry alone or in combination with first coordination sphere descriptors. Apart from demonstrating metal sensitivity and yielding better benchmarking results, this new version allows the assessment of the influence of considering the metals first coordination sphere versus backbone pre-organization on how metallic species bind to proteins.

Carotenoids are structurally diverse isoprenoid pigments that play central roles in photosynthesis, photoprotection, membrane organization, and cellular signaling. Despite their biological and technological importance, atomistic simulations of carotenoids remain limited by the lack of a transferable force field spanning the chemical diversity of naturally occurring compounds, including glycosylated and acylated derivatives. Here we present CARS (CARotenoidS), a transferable force field for carotenoids that integrates seamlessly with the AMBER family of biomolecular force fields. Parameters were systematically optimized against 22 957 r2SCAN-3c reference energies for 25 representative molecular fragments, yielding an accurate description of polyene conformational energetics, ring rotations, and molecular geometries. Across diverse validation sets, CARS substantially outperforms GAFF2 and provides improved agreement with quantum-chemical reference data for glycosylated and acylated carotenoids. For zeaxanthin, CARS also surpasses OPLS-AA, CGenFF, and previously published carotenoid-specific parameters set in reproducing conformational energetics and structural properties. Two complementary parameter sets are provided: CARS for glycosylated and non-lipidated carotenoids, and CARS+Lipid21 for carotenoids containing saturated or monounsaturated lipid chains. By providing the first unified and transferable parameterization covering the structural diversity of natural carotenoids while remaining fully compatible with established AMBER force fields, CARS removes the need for molecule-specific reparameterization and enables reliable molecular simulations of carotenoids in proteins, membranes, and other complex biological assemblies.

Numerous advances in quantum and computational chemistry over the last decades, well as the development of computer science, allowed utilisation of more precise and complex models, which can be now applied to much bigger systems than in the past. The authors used Gaussian, coupled with theoretical methods, to predict a new way of peptide bond formation, which could have taken place in prebiotic conditions. To better tackle this difficult task, the properties of substrates (glycine-derived radicals) were extensively analysed, using the aforementioned tool - Gaussian, paired with taking resonance and hybridisation into account, to better understand the stereochemistry and the very nature of processes taking place. The result is a series of reactions, which without any sophisticated catalysts and with relatively low energy thresholds ({inverted exclamation}20 kcal/mol) can lead to formation of dipeptides (and further, oligopeptides). The authors also hope, the other predicted properties of the investigated molecules can be of use to any researcher, who would like to utilise them in their experiments. Author summaryOur goal was to investigate a way first peptide bonds in prebiotic conditions could have been formed. This is an extremely important step in research into the beginning of life on Earth. We found a very promising series of reactions, which uses atomic hydrogen as its only catalyst and confirmed our expectations with theoretical calculations, using Gaussian. There are two radicals derived from glycine, which perform major roles in the process, so we investigated their properties with Gaussian and verified that the results are in agreement with our own theoretical considerations. This involved checking for possible geometric isomers and conformers and creating models which could explain their properties. We are well aware that such calculations have limitations and there is no model, which is 100% accurate, so our results should be further confirmed by empirical data in the future. However, we still to be as thorough as possible in how we approached the subject.

The analysis of nucleic acid secondary structures is overwhelmingly dominated by methods that analyze the thermodynamic equilibrium distribution and which ignore all dynamic aspects of nucleic acid folding. Yet, there are numerous popular examples of nucleic acid folding that rely on kinetic models, such as RNA riboswitches or DNA strand displacement systems. Here, I am presenting fuzzyfold, a Rust-based software package for nucleic acid secondary structure analysis with an explicit focus on stochastic modeling. The framework introduces three-way and four-way shift moves with a biophysically motivated rate-model parameterization, and it is developed with an emphasis on both model flexibility and performance, e.g. allowing for the generation of single co-transcriptional trajectories for thousand-nucleotide long RNA molecules in just a few minutes. The main strength of the fuzzyfold package, however, is its focus on user and developer interfaces for long-term development. It provides easily installable command-line interfaces, e.g. for aggregating data from multiple parallel trajectories efficiently into an ensemble-level dynamic analysis. For developers, the code-base supports straight-forward substitution of thermodynamic and kinetic free-energy models, and a flexible library interface with Python bindings, enabling integration of individual components into custom computational workflows.

Adaptive sampling accelerates the exploration of conformational space in molecular dynamics (MD) simulations by repeatedly analyzing the accumulated trajectories and seeding a new round of simulations from informative configurations. A growing collection of adaptive sampling policies has been proposed, each built around a particular notion of what makes a configuration informative, yet these methods are scattered across separate and often incompatible implementations, which complicates their systematic comparison and their combined use in meta adaptive sampling schemes. Here, we present AdaptivePy, a compact and extensible Python framework that implements nine seed-selection policies behind a single configuration-driven interface, spanning simple population-based baselines, several established machine-learning and geometry-based methods, and two ensemble or meta sampling policies introduced in this work. We show that the shared implementation reproduces the characteristic selection behavior of each policy on a series of analytic benchmark landscapes. We also introduce a new adaptive sampling scheme that employs TS-DAR, a deep learning framework originally designed to identify transition states, into an acquisition criterion that drives the discovery of an entire multi-basin landscape starting from a single basin. We further demonstrate that the common interface enables meta adaptive sampling policies, which aggregate the rankings of several policies into a single set of seeds. AdaptivePy thereby provides a unified testbed for the adoption, benchmarking, and continued development of adaptive sampling methods for biomolecular MD simulations.

Aptamers are structured nucleic acid ligands capable of high affinity, high specificity molecular recognition generated using variations of the SELEX (Systematic Evolution of Ligands by Exponential Enrichment) process. However, SELEX often produces sequences that enrich yet may lack binding efficacy. We propose a measure called the Ruggedness Composite Index (RCI) along with a method for computing it, that can be used to distinguish binding-competent ('active') aptamers from weak or non-binding ('inactive') aptamers. Given a set of aptamers, RCI incorporates information on their fragmentation (landscape partitioning), basin entropy (metastable state distribution), cumulative density irregularity (non-uniform occupancy), and structural energy correlation length (structure-energy coupling scale). We test whether secondary-structure folding energy landscape topology distinguishes active from inactive aptamers using a multiscale level set framework across six datasets. Active aptamers show lower RCI values and occupy smoother, funnel-like conformational spaces, while inactive aptamers show higher RCI values, reflecting fragmented, high-entropy landscapes. By contrast, classical thermodynamic features, such as minimum free energy, show limited discrimination between active and inactive aptamers. In all datasets, sequences that exhibit enrichment which is not monotonic but lack specificity exhibit elevated ruggedness, indicating landscape topology can predict non-specific enrichment. These results indicate that folding landscape organization can be used as a predictor of aptamer activity and establish RCI as a simple, mechanistically interpretable measure for improving candidate prioritization, especially in therapeutic aptamer discovery.

The molecular mechanics-generalized Born surface area method (MMGBSA) is one of the most commonly used end state approaches used for the calculation of the binding free energy towards computational drug design and screening studies. It is customary to break up the free energy into van der Waals, electrostatic, polar solvation (GB), and nonpolar solvation (SA) terms and then either correlate these terms with experiment or assign physical meaning to each term. Here, we demonstrate that this assumption of independent fitting coefficients for decomposed energy terms could be invalid. Through analytic derivation and large-scale molecular dynamics simulations, we show that (i) the protein and ligand Coulomb interaction energy and the GB solvation correction are almost perfectly collinear (R2[≥]0.99) reflecting their designed role as vacuum electrostatics plus solvent screening, and (ii) the van der Waals interaction and SA term likewise exhibit strong correlation, as both depend primarily on buried surface area. Interaction entropy and C2 entropy corrections are also found to be strongly dependent on underlying electrostatic fluctuations, further reinforcing redundancy. These findings hold both at the level of instantaneous trajectory fluctuations and when averaged across a diverse set of 139 protein-protein complexes and persist in both single-trajectory and three trajectory MMGBSA protocols. Our results caution against using decomposed MMGBSA terms as independent predictors in regression models and suggest instead combining correlated terms into effective polar, nonpolar, and entropic contributions. Our study provides a systematic diagnosis of collinearity in MMGBSA and highlights pathways toward more interpretable and statistically robust predictive modeling.

We used simulations of short single stranded DNA (ssDNA) homopolymers, based on the sequence dependent Three Interaction Site (TIS) model, to calculate the persistence length (lp) in polyvalent cations. The TIS model accounts for stacking interactions and electrostatic interactions are treated using the Coulomb potential. We find that lp for dT30 (T is thymine) and dA30 (A is adenine) is quantitatively fit using [Formula] ([Formula] is the bare persistence length,{lambda} is a dimensionless constant, and{kappa} -1 is the inverse Debye length) in the divalent cations Mg2+ and Ca2+. The dependence of lp on{kappa} is surprising because it was derived for long flexible polyelectrolytes in which the charges interact via the Debye-Huckel potential. The [Formula] values are 0.4 nm and 1.1 nm for polyT and polyA, respectively. Strikingly, lp is almost independent of the tetravalent spermine concentration. There is no clear theoretical explanation although simulations suggest that the number of spermine molecules that bind to the ssDNA saturates at a small value. A qualitative picture, based on the restrictions of access to the phosphate groups due to volume exclusion of the anisotropic structure of Spm4+, rationalizes the simulation results. The predicted dependence of lp in spermine awaits experimental test.

Interleukin-13 (IL-13) is an immunomodulatory cell signaling cytokine that has been implicated in neurodegenerative disease and chronic inflammation. IL-13 binds to its low and high affinity receptors, IL-13 receptor 1 (IL-13R1) and IL-13 receptor 2 (IL-13R2), respectively, with residence times that vary accordingly. As the binding kinetics of the cytokine-receptor complexes influence cellular responses, we employed the molecular dynamics (MD) simulation-based{tau} -random acceleration molecular dynamics method ({tau}RAMD) to compute relative residence times for wild-type (WT) IL-13 and 19 IL-13 mutants to the two receptors. Comparison with experimental kinetic data shows that the{tau} RAMD computations capture the trends in residence times. Analysis of simulated dissociation trajectories of the cytokine-receptor complexes reveals two distinct dissociation pathways of IL-13 from each of the receptors. This study thus pinpoints key determinants of the interaction of IL-13 with its receptors which could be targeted for therapeutic design. Statement of SignificanceCytokines are regulatory proteins that bind to cell surface receptors and thereby send signals to the cellular interior. Interleukin-13 (IL-13) is a cytokine that has a low and a high affinity receptor. It has important physiological roles, and its deregulation is involved in diseases such as atopic dermatitis and asthma. We employed a molecular dynamics simulation-based method to compute the effects of changes in the sequence of IL-13 on the lifetimes of complexes of IL-13 and its receptors. Comparison with experiments supports the validity of the computational approach and analysis of the simulations reveals two distinct ways in which IL-13 dissociates from each receptor. These results thus provide a map for targeting IL-13 - receptor interactions for the design of therapeutics.

The process of drug unbinding is of immense importance in the field of biophysics and therapeutics. The behavior of these systems is greatly influenced by their thermodynamic and kinetic properties. Therefore, it is crucial to accurately estimate the ligand binding free energies and rate of ligand dissociation, yet these processes are often governed by rare event transitions that lie beyond the reach of standard brute-force molecular dynamics simulations. While enhanced sampling simulations offer a solution, their efficacy is strictly contingent upon the selection of appropriate collective variables (CVs) which is non-trivial for complex systems like protein-ligand complexes. In this study, we present a method to derive optimized CV from transition state region (TS) via an interpretable machine learning (ML) model, Elastic Net. By employing some physically intuitive order parameters, the derived optimized CV from the TS-region greatly accelerate ligand binding-unbinding transitions and achieves rapid free energy surface (FES) convergence across diverse systems including buried and solvent exposed active sites such as Trpsin-benzamidine complex, host-guest systems and sodium epoxidase etc. Intriguingly significant contribution of the ligand hydration is found in the optimized CV which depicts crucial role of solvent in driving ligand binding-unbinding transitions. The estimated binding free energies for different protein-ligand complexes match quite well with experiments, while maintaining a low computational cost. The derived optimized CV is also used to calculate the ligand residence times across different systems and calculated residence times are within the experimental range for all systems, again with very little computational costs. Moreover, we show that the optimized CV constructed from TS region via an interpretable ML model is transferable across diverse systems, offering a robust and scalable framework for drug discovery and investigation of complex biomolecular recognition.

Alzheimer's disease is characterized by the accumulation and aggregation of amyloid-{beta}(A{beta}), but the molecular mechanisms linking environmental and infectious factors to A$\beta$ conformational changes remain incompletely understood. Herpes simplex virus type 1 (HSV-1) has been proposed as a potential contributor to AD pathology, and interactions between the viral glycoprotein B (gB) and A$\beta$ may influence the conformational behaviour of the peptide. Molecular dynamics (MD) simulations provide atomic-scale information on such interactions, but conventional structural descriptors may not fully capture changes in the organization of residue interaction networks. Here, we introduce a graph-geometric framework based on Forman-Ricci curvature to characterize the evolution of residue interaction networks during MD simulations. Each simulation frame is represented as a residue interaction graph based on C--C contacts, and residue-wise curvature profiles are analysed across time. We apply the framework to A{beta}1-42 in isolation and in complex with HSV-1 gB. Conventional MD analyses indicate stable association of the simulated complex, favourable interaction energetics, and conformational changes in A{beta}, including a transition from -helical structure toward {beta}-turn-rich conformations over the simulated timescale. Forman-Ricci curvature reveals pronounced and spatially localized remodelling of the A{beta} residue interaction network in the complex, with the strongest changes concentrated in the C-terminal region. These regions also exhibit reduced temporal curvature fluctuations and progressively distinct geometric behaviour throughout the simulation. Hierarchical clustering further identifies cooperative groups of residues with coordinated curvature dynamics, including a prominent C-terminal domain. Together, these results demonstrate that Forman-Ricci curvature provides a complementary description of biomolecular dynamics by capturing changes in the geometric organization of residue interaction networks that are not directly represented by conventional structural descriptors. The framework provides a general computational approach for studying network-level structural remodelling in protein molecular dynamics and offers a quantitative perspective on the conformational consequences of HSV-1 gB--A{beta} association.

Neutron scattering is a powerful technique to investigate atomic structures and molecular dynamics of proteins at the nano-scale. When it comes to dynamics, incoherent and coherent scattering respectively provide information on the single and collective dynamics of nuclei. In proteins, hydrogen has the highest incoherent cross-section, and it is common practice to overlook the contribution of coherent terms stemming from all nuclei. However, the fast collective dynamics of heavier nuclei could also be studied if coherent scattering and incoherent scattering were experimentally separated. The recent advent of polarized neutron spectroscopy with sufficient flux and energy resolution has made it possible, and opens new perspectives to investigate the relative importance of coherent scattering and the information it provides on biological samples. The present study reports on the use of polarized quasi-elastic neutron scattering (QENS) and the application of a minimalistic model adapted to both individual and collective dynamics. Using a perdeuterated green fluorescent protein as a model globular protein, the study provides an interpretation of the dynamical parameters obtained with QENS, and a comparative study of the Elastic Coherent and Incoherent Scattering Factor. Based on both experiments and calculations, we discuss the relative importance of distinct and self components of coherent scattering, which is often wrongly assumed to be representative of collective dynamics only. The results highlight the current impediments rendering complicated a straightforward analysis of fast collective dynamics in hydrated protein samples.

Biomolecular condensates formed by intrinsically disordered proteins require molecular models that accurately describe proteins in both dilute solution and condensed phases. Explicit-solvent coarse-grained models offer an attractive balance between chemical resolution and computational efficiency. Yet, it remains unclear whether improving dilute-state properties is sufficient to obtain an accurate description of condensates. Here, we address this question by introducing minimal modifications to the Martini 3 force field that combine recent advances in bonded interactions with refined protein-water interactions and strengthened glycine self-interactions, while preserving the underlying chemical transferability of the model. The resulting model substantially improves the description of single-chain conformations across a diverse benchmark of disordered proteins. We then investigate phase separation of the well-characterized low-complexity domain of heterogeneous nuclear ribonucleoprotein A1 and its sequence variants. The model reproduces several key physicochemical properties of biomolecular condensates, including chain expansion in the dense phase, sequence-dependent intermolecular contacts, protein diffusion and its relation to single-chain dimensions, and hydration, while revealing quantitative limitations in condensate density, phase equilibria, and ion partitioning. Our results show that improving dilute-state behaviour translates into a better description of condensed-phase properties, including condensate density, but is not sufficient to quantitatively reproduce the equilibrium between the dilute and dense phases.

Background: RNA molecules form thermodynamic ensembles, but interpretation often requires a single representative structure. Existing base-pair centroid estimators assess agreement at the level of individual base pairs and do not directly target nesting depth along the sequence. Methods: We introduce Mountain Centroid, which minimizes expected squared mountain-profile distance, and derive dynamic programming algorithms with and without RNA pairing constraints. We also combine the Mountain Centroid objective with the base-pair centroid gain. Results: Across 21,254 RNAStrAlign sequences, Mountain Centroid had lower median normalized mean squared mountain distance (NMSMD) than minimum-free-energy (MFE) and base-pair centroid ({gamma} = 1) structures, whereas its median base-pair F1 was lower. Imposing RNA pairing constraints improved base-pair F1 for 59.35% of sequences and reduced it for 3.58%. At an illustrative weight, the combined objective had median base-pair F1 similar to MFE while retaining lower median NMSMD than MFE and all tested {gamma}-centroid settings. Conclusions: Mountain Centroid represents an RNA structural ensemble with a single secondary structure that reflects how nesting depth varies across nucleotide positions. Combining mountain-profile and individual-base-pair criteria allows their relative contributions to be varied.

Martini 3 is a force field ideally suited to simulating long intrinsically disordered proteins (IDPs) in cell-like surroundings. So far, most Martini 3 variations intended for IDPs have only been benchmarked on shorter IDPs of up to 140 amino acids. In this paper, we present a comprehensive benchmark including IDPs up to 809 amino acids in length and compare the behavior of four well-known Martini 3 variations for IDPs. Modifications to only the bonded parameters result in excessively compact conformations, thereby failing to reproduce the experimental radius of gyration observed for large IDPs. In contrast, general rescaling of interaction parameters, including tuning electrostatic interactions in the case of highly-charged long IDPs, yields acceptable levels of compaction at all tested length scales.

Quantifying how cooperative, many-body relationships drive allostery in protein networks remains a major challenge. To address this, we develop the Laplacian minor hierarchy, a mathematical framework that characterizes the geometric invariants of a protein network. Lower-order minors yield standard metrics including the partition function and effective distances, whereas higher-order minors define novel topological measures: cooperation indices, each bounded between zero and one, that characterize pathway correlations at increasing levels of complexity, the third-order minor determines whether allosteric pathways are correlated or uncorrelated, and the fourth-order minor quantifies how distinct pathways communicate through intermediary residues. We apply this framework to analyze the evolutionary adaptation of the PSD95pdz3 domain from Class I to Class II ligand specificity via mutations G330T and H372A. The cooperation index demonstrates a distinct evolutionary hierarchy: the G330T mutation establishes distributed pathway couplings that the H372A mutation subsequently exploits, whereas H372A alone produces minimal global changes. Furthermore, the fourth-order analysis identifies His317 as a critical intermediary node bridging the class-switching (330-372) and class-bridging (330-400) allosteric pathways. These results demonstrate that allosteric dependencies emerge only when mutations accumulate in specific combinations, with a hierarchical organization of pathways structured around position 330 and intermediary nodes His317 and Phe400. Rather than predicting allosteric mechanisms, this framework provides a mechanistic explanation for why and how allostery emerges during protein evolution.

Accurate prediction of Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) properties is a central challenge in early-stage drug discovery, where experimental determination remains costly and time-consuming. In this work, we propose a quantum-inspired preprocessing framework in which statistical dependencies among molecular descriptors are encoded into a parameterised many-body Hamiltonian, and the expectation values obtained by simulating its time evolution serve as additional inputs to a gradient-boosted ensemble model (CatBoost). Mutual information (MI) is used both to select the most informative descriptors and to set the coupling strengths of the Hamiltonian, so that the induced entanglement structure reflects empirically measured feature correlations; the evolution is realised with a short digitised-counterdiabatic schedule that generates a compact set of expectation-value features while keeping the circuit shallow. The resulting quantum-derived feature vectors are concatenated with the full MapLight descriptor set, concatenated ECFP, Avalon, and ErG fingerprints together with RDKit physicochemical properties, before training. We evaluate the pipeline on the AqSolDB aqueous solubility benchmark from the Therapeutics Data Commons (TDC) platform, achieving a mean absolute error (MAE) of 0.746 +/- 0.006 log(mol/L), which is within the reported error bars of the current top-performing model on the TDC leaderboard (MAE = 0.741 +/- 0.013). Ablation experiments show that the quantum-derived features match classical second-degree polynomial interaction features derived from the same MI-selected subset, while forming a far more compact representation (85 quantum features versus up to 4,950 polynomial terms, an approximately 58-fold reduction). SHapley Additive exPlanations (SHAP) analysis identifies the physicochemical drivers of solubility predictions, offering interpretable insight into model behaviour. These results demonstrate that MI-guided Hamiltonian feature extraction can reproduce the performance of strong classical interaction models on aqueous solubility while generating a compact, interpretable feature representation that is compatible with future quantum execution.

Many functional RNA motifs adopt structures that deviate from the canonical A-form helix and are emerging targets for RNA-directed therapeutics. The microtubule-associated protein tau (MAPT) A-bulge motif (5'-GCAGU/5'-ACGU) is one such motif. Because its structure is stabilized by a delicate balance of local interactions, its accurate modeling remains a major challenge for molecular dynamics (MD) simulations. The experimentally determined nuclear magnetic resonance (NMR) structure of the MAPT A-bulge motif provides a stringent test of whether RNA force fields can accurately reproduce the experimentally observed conformation. Most current AMBER-family RNA force-field models have incorrectly favored a non-native base-triple state of the MAPT A-bulge motif over the experimentally observed stacked state. Structural comparison of the stacked and base-triple conformations revealed that overly favorable NH-N hydrogen bonds between the bulged adenosine and an adjacent Watson-Crick base pair were the primary source of this imbalance. We developed gHBfix-18Ab, an 18-component hydrogen-bond correction that distinguishes NH and NH2; donors. gHBfix-18Ab was combined with the previously developed OL3CP and NBfix0BPh corrections to generate the composite model gHBfix-18Ab*. This model restored the experimentally observed stacked state as the global minimum in the calculated free-energy profile and improved agreement with NMR-derived distance data for the A-bulge region. Importantly, gHBfix-18Ab* did not produce marked structural destabilization of the cUUCGg tetraloop, a widely used benchmark for RNA force-field validation, suggesting that the refinement preserves the stability of the unrelated RNA motif. These results demonstrate that targeted refinement of hydrogen-bond interactions provides a practical strategy for systematic improvement of RNA force fields toward more accurate modeling of noncanonical RNA motifs.

Molecular dynamics simulations of nucleic acids are performed using a solvent-buffer distance of 10 [A] between the solute surface and the simulation box boundary. Although this cell size has been extensively explored in protein simulations, its implications for nucleic acid dynamics are not well understood. Nucleic acids are elongated, highly charged, and flexible structures with hydration and dynamical properties distinct from those of proteins and therefore, they may require different solvent-layer considerations in simulations. In this study, we investigated the effect of simulation cell size on nucleic acid dynamics by simulating a 30-base-pair double-helical nucleic acid structure and its two single-stranded forms using solvent-buffer distances of 3, 5, 10, 15, and 20 [A]. Smaller cells may impose restricted hydration, molecular crowding, and periodic image interactions. However, larger cells provide solvent space for conformational relaxation. A total of 45 s of molecular dynamics simulations were performed (3 structures x 5 cell sizes x 3 replicates x 1 s). Our results show that while the commonly used 10 [A] buffer may be sufficient to maintain the stability of the double-stranded nucleic acid, larger cells are required to capture the conformational dynamics of single-stranded structures. In both, increasing the cell size to 15 or 20 [A] enables broader conformational sampling. The first hydration shell exhibits reduced crowding in the 20 [A] cell, consistent with more relaxed conformations. At larger cell sizes, single-stranded nucleic acids adopt compact, self-associated conformations for stability. Together, this study presents physical insight into how simulation cell size and solvent environment influence nucleic acid dynamics.

Circuit Topology (CT) provides a fundamental framework for analysing folded polymer chains, with applications in functional annotation, protein engineering and drug development. We present a protein CT analysis plugin for PyMOL v3.1.6.1 with a graphical user interface (GUI), automatic installation, and novel features developed through integration with PyMOL's application programming interface (API). The plugin integrates various previously developed CT methodologies for studying structured proteins and their complexes as well as the dynamics of disordered proteins. Analysis of a representative protein and a molecular dynamics trajectory demonstrates the plugin's three analysis modes and their outputs. The plugin reproduces the reference ProteinCT implementation exactly on the structures tested, and is distributed with a versioned release, a pinned environment and a one-command reproduction of every result reported here.